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Abstract. We present a purely algebraic formulation (i.e. polynomial equations only) of the minimum- 
cost multi-impulse orbit transfer problem without time constraints, while keeping all the variables with 
a precise physical meaning. We apply general algebraic techniques to solve these equations (resultants, 
Grobner bases, etc.) in several situations of practical interest of different degrees of generality. For 
instance, we provide a proof of the optimality of the Hohmann transfer for the minimum fuel 2-impulse 
circular to circular orbit transfer problem, and we provide a general formula for the optimal 2-impulse 
in-plane transfer between two rotated elliptical orbits under a mild symmetry assumption on the two 
points where the impulses are applied (which we conjecture that can be removed). 


1. Introduction 

Since the start of the space age with the first launch of a satellite to space (Sputnik in 1957), the interest 
in the study of space maneuvers that use the available resources like time and fuel efficiently has been 
growing steadily. In many real life cases, a satellite can serve multiple purposes, requiring for that a 
change of its orbit. Some maneuvering is also needed during the initial launch of a satellite, or when a 
spare satellite has to be brought to its intended orbit. 

Maneuvering a satellite can be done in two different ways: continuous thrust or a sequence of instanta¬ 
neous and discrete impulses. This paper focuses on the latter. 

The orbit transfer problem with a fixed time of flight was studied by Lambert, who provided a solution 
in the case of two impulses. For a discussion of this problem, see m- 

However, the scarcest resource in space is fuel, since it represents a load on the spacecraft that cannot 
be too large to avoid launching problems and to reduce costs. For this reason, we focus our attention on 
the minimum fuel transfer problem with unconstrained time. Using the well-known Tsiolkovsky rocket 
equation, we consider the sum of the individual impulses (difference between velocities before and after 
the thrust is applied) as the cost function. This usually appears in the literature as Av. 

In the case of a transfer between two circular coplanar orbits, Hohmann gave an explicit solution with 
two impulses in [7], which was later proven optimal analytically by Barrar [3]. For the case with three 
impulses, Hoelker and Silber [5] have shown that a bi-elliptical transfer has a lower fuel requirement 
than the Hohmann transfer for some special initial and final orbits. Roth m extended the notion of 
bi-elliptical transfer to the case of two inclined orbits. 

In this paper, we provide a detailed study of transfers between two circular orbits, including out-of- 
plane maneuvers and also the possibility that the initial and final angular momentum point in opposite 
directions. In all cases, we have proven algebraically that the Hohmann transfer is optimal for two 
impulses. 

Another problem of interest is transferring a satellite between predetermined points in the initial and 
final orbits. This situation was studied by Avendaho and Mortari in where they provided a closed- 
form solution. Here, we reobtain this solution algebraically, applying a much more efficient method. 
Previous attempts to solve this problem involved the use of iterative methods or an equation that needs 
to be solved numerically (see [21 [SI ISl HI US] ) ■ 
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The last problem we study is the optimal transfer between two identical ellipses that are coplanar and 
rotated a fixed angle. This case was studied numerically by Bender in [3]. However, we provide an 
algebraic solution, which is fully explicit under a mild symmetry assumption. 

This paper is organized as follows. In SectionJ^we have compiled all the equations of Celestial Mechanics 
that we will need in the paper. In Section we present an algebraic approach to the multi-impulse 
minimum-cost orbit transfer problem. We have put special emphasis in explaining the physical meaning 
of all the variables involved. Three kinds of problems are studied: point to point, point to orbit and 
orbit to orbit. We also consider two possible cost functions. 

In Section we present our solution to the general point-to-point problem with two impulses and cost 
function as in In Section]^ we provide a solution to a generalized version of the Hohmann transfer 
where out-of-plane maneuvers are allowed. 

Finally, the orbit-to-orbit problem between identical and coplanar orbits is studied in Section]^ The 
solution we obtain requires solving a large system of polynomial equations, which can be solved explicitly 
if we assume a symmetry condition. We have done extensive numerical tests showing that the symmetry 
condition is always satisfied in them. 


2. KePLERIAN MOTION 


The motion of a particle in a Keplerian gravitational force field is given by the solution of the second- 
order differential equation 


( 1 ) 


r{to) = To, r{to) = uo, r = 


where /r > 0 is the standard gravitational parameter of the field, tq and Vq are the initial position and 
velocity, and r{t) is the position of the particle as a function of time. For any solution of Eq. 0, the 
angular momentum vector 

h = fx r, 

the eccentricity vector 

( 2 ) - " 


e = 


M 1^1’ 
Iri^ 


and the total energy 

I ^ O 

A* 

" 2 |rl 

are constants with respect of time El Ch. 8.3]. The vectors h and e are always orthogonal, i.e. 

h - e = 0, 

and any pair of mutually orthogonal vectors h and e can be obtained for some initial conditions tq and 
uq. Moreover, the total energy satisfies 


I-|ef = - 


2E\h\'^ 


so its value can be determined from h and e alone when h ^ 0. 
The angular momentum is always orthogonal to r, i.e. 

(3) r-h = 0, 

so the motion is planar. Besides, it follows from Eq.0 that 


( 4 ) 


^ fx h _ j/iP 

\f\ + e-r= -r= - 

/i n 
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is also constant, which is the implicit equation of a conic with one focus at the origin, eccentricity e = |e| 
and semilatus rectum p = |/ip//r, when p 7 ^ 0. In the case e > 1, i.e. when the conic is a hyperbola, the 
equation describes only the branch in which the particle is moving. The degenerate case p = 0 will be 
discussed later in this section. Finally, multiplying Eq. (§ by h and moving some terms, we obtain an 
expression for the velocity of the particle at any given position: 


( 5 ) 


r = 




1 ^ 


It is important to note that, when the orbit is an ellipse (e < 1), any point r that satisfies Eq. ([^ and 
Eq. 0 will be visited by the particle at some time t > to since the motion is periodic. However, this is 
not true when e > 1 because in the case of a parabolic (e = I) or hyperbolic (e > 1) trajectory only the 
points satisfying the extra condition 

(f X fo) ■ h < 0 


will be visited. 

As we mentioned above, the case p — 0 needs to be discussed separately. Here we have h = 0, so rg 
is parallel to vq- The eccentricity vector e = — ypy is constant, so the trajectory is contained in the 
line through the origin with direction e. There are two possible cases: either the initial velocity is high 
enough to escape the gravitational attraction of the held or the particle will hrst move in the direction 
Vo until a point where its velocity becomes zero and then come back towards the origin, thus entering 
in a periodic motion. 

To avoid the extra complexity needed to handle parabolic and hyperbolic motions, as well as the degen¬ 
erate case p — 0 described above, we restrict our analysis to elliptic orbits, e < 1 , and non-degenerate 
trajectories, h ^0. 

In order to work with polynomial equations, we need to remove the divisions, the square roots and the 
constant p from some of the equations above, so we introduce the vectors 


% r f p h ^ p _ 

r =— IT, w = I = s = i X e. 

1^ Vm \h\^ 


Note that I and s are orthogonal, i.e. 

( 6 ) 


I • s = 0, 


and that the angular momentum and eccentricity vectors can be simply recovered as 

r ^ I ^ s X I 

h = Jp—, e = 

'^ 1712 ’ 1712 


Of course, the case I — 0 has to be excluded, and the condition e < 1 translates into [sj < |Z|. 

Any unit vector r orthogonal to I, i.e. 

(7) r • r = 1, f-1 = 0, 

determines a point on the orbit. The exact location can be obtained from Eq. 

( 8 ) ^ = + 

and the velocity of the particle at that point is, according to Eq. 0, 

(9) ut = s + I X f. 

Einally, note that in the case of elliptic orbits (|^ < |^^), the right-hand side of Eq. 0 is always positive, 
so no extra inequalities are needed to guarantee a valid value of 
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Figure 1. Ellipse with focus F and a satellite S on it. 


3. Multi-impulse orbit transfers 


An n-impulse orbit transfer is represented algebraically by the vectors 

/q 5 ^ 1, ■ • ■ ; : "^0; ^ 1J • ■ ■ J 5 '^0 ; j • • ■ ; ^n— 1 

where the pair {li,Si) determines the i-th orbit and ry corresponds to the point where the impulse is 
applied to change from the i-th to the {i + l)-th orbit. These vectors, according to Eqs. (|§,0 and 
are constrained by 

( 10 ) k- Si = 0 , 

( 11 ) k ^ 0 , 

( 12 ) 

for z = 0 ,..., n, and 


I Si I ^ l^i IJ 


(13) 

(14) 

(15) 

(16) 


li ‘ Ti — 0 , 

^i-t -1 ‘ — 0 , 

= 1 , 

l^i| ~b (Si ^ ^i) ‘ I’ll — “t (Si-t-1 ^ ^i-t-l) ' 


for alH = 0,..., n — 1. Conversely, any sequence of vectors satisfying all these restrictions represents a 
valid n-impulse transfer. Moreover, all the equations are invariant under rotation by a fixed angle and 
rescaling of the li and 

The following table shows the number of unknowns and algebraic equations that dehne the configuration 
space for each type of n-impulse transfer (n > 2 ). 



3d-transfer 

2 d-transfer 


^unknowns 

^equations 

^unknowns 

^equations 

Point to Point 

9n- 12 

5n- 5 

5n- 7 

2n — 2 

Point to Orbit 

9n — 9 

5n- 3 

5n- 5 

2n- 1 

Orbit to Orbit 

9n — 6 

5n- 1 

5n — 3 

2n 


In the orbit to orbit problem, the vectors Iq, Sq, In, Sn are given and the remaining variables are 
considered unknowns. In the point to orbit, the initial point is given, so Tq is also known, thus reducing 
the number of unknowns (and equations). Finally, in the point to point problem, the final point is also 
given, i.e. r„_i is known. 
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The two-dimensional version of these problems considers all orbits in the z = 0 plane, so h = (0,0, kz), 
Si = (six,Siy,0) and fl = (xi,yi,0), hence the reduced number of variables and equations needed to 
handle them. 

The velocities at the points fl, immediately before and after the impulse is applied, are written w, and 
w *, respectively. It follows from Eq. ([^ that 

(17) Wi = Si + li X ri, 

(18) w* = Si+i + k+i X fi. 


The cost (fuel-wise) of such a transfer is proportional to the sum of = Iwi — w;*|, denoted hereafter 
by fi- To avoid the square roots that are implicitly present in Aj, we also consider a cost function /2 
which is the sum of the squares of the A^: 


-»* 12 
w. 


n—1 n—1 

fl = '^\Wi-W*\, f2 = '^\Wi-‘ 
i=0 i=0 

If the vectors li and s) are rescaled by a factor c, then fi and /2 are multiplied by a factor |c| and |cp, 
respectively. 

When the cost function /i is used, the trick to avoid the square roots consists of considering A^ as a 
variable, redefining the cost function as 

n—1 






and adding the algebraic equations 

~ ^i+ll X (^li , 

for z = 0,... ,n — 1. The last equations can be obtained by substituting Eqs. 0 and ( |I^ into the 
definition of A^ : 


A^ =\w,-w*\ = |si + liXn - s,+i - I 


i+l 


X r. 


— I \^i ^I+l| T ^((si Si-|-i) X (^li ^^_|_i)) 


At this point we have a classical problem of constrained minimization, which we approach with Lagrange 
multipliers. 


Theorem 1 (Lagrange multipliers). Let q,qi,... ,qm ■ K* —t R in C°° and p G a common zero 
of qi,... ,qm be such that the vectors Vgi(p),..., V(7 m(p) are linearly independent. Then p is a local 
extremum of q on the manifold defined by {qi = • • • = 9 m = 0} if and only if there exists Ai,..., Am G K 
such that "S/q(p) = AiV 9 i(p) \rn^<lm{p)■ 


In our case, we have an algebraic variety V = {gi = ... = 9 m = 0} C K^, defined by polynomials 
9 i,..., 9m G M[xi,..., Xfc], and another polynomial function q that we want to minimize on V. In order 
to apply Theorem]^ we need to exclude first the points where Vqi{p),... ,Vqm{p) are not linearly 
independent, which is, by definition, the set of singular points V* of V. Computationally, V* is the set 
of points of V where all m x to minors of the matrix [dqi/dxj]i<i<m,i<j<k have zero determinant: 


W= pGM'^ : 9i = ---=9m = 0 A ''' ’ = 0, VJ = {ji, ■ • •, Jm} C {I,..., fc} 

I d{xp,...,XjJ 

These points have to be considered critical points (i.e. they are potential local extrema) and have to be 
evaluated separately. 

On the remaining points, E \ E*, the local extrema can be found directly by Theorem]^ solving the 
system of to -|- fc equations 91 = • • • = 9 m = 0 and V 9 = A 1 V 91 -I- • • • -I- AmVgm in the m + k unknowns 
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a;i,...,Xfc, Ai,..., Am G K, and disregarding the solutions with (xi,... ,Xk) G V*. Removing these 
solutions is actually not needed since they are always critical points. The set of solutions of the m + k 
equations described above and the set of all critical points of q are denoted Vq and respectively: 

Vq = {(xi,..., Xfc, Ai,. ■ •) Am) G : qi = ■ ■ ■ = qm = 0 A Vq = AiVgi + • • • + AmV^m} ) 

(19) = R* U (14) C 

where tt/j : —>• is the projection onto the first k coordinates. 

Including the Lagrange multipliers, and the extra variables for i = 0,..., n — 1 when minimizing fi 
instead of / 2 , we get the following total number of unknowns (which is equal to the number of equations): 



min 

3d-transfer 

(/i) 

2d-transfer 

min 

3d-transfer 

(/2) 

2d-transfer 

Point to Point 

Point to Orbit 
Orbit to Orbit 

16n - 17 
16n - 12 

16n — 7 

9n — 9 

9n — 6 

9n — 3 

14n - 17 
14n - 12 
14n- 7 

7n — 9 

7n — 6 
7n-3 


By using standard linear algebra, it is possible to eliminate all the Lagrange multipliers: 


( 20 ) 


xrcrit _ xr* 

Vq — V 


u 


qi = ■■■ = qm = 0 
d{q,qi,...,qm) 




= 0, VJ = {ji,...,jm-Li} C {!,..., A:} 


The expression above shows that Vq^'^* can be written as the union of two algebraic varieties in . 


4. Minimum Av'^ Lambert problem 


In this problem, the vectors rg, rq, wq, wl are known, from which the vectors Iq, I 2 , soj S 2 can be 
computed directly. The unknowns are li and ^ 1 , from which we can deduce and wi. There are two 
cases, depending on whether rg and rq are linearly independent or not. 

In the first case, we can assume without loss of generality that rg and rq both lie on the ccy-plane, so the 
unknowns can be written li = (0,0, liz) and ?l = (sia:, siyj 0). We can further assume that rg = (1,0, 0) 
and rq = (a;i,i/i,0) with yi ^ 0 and x\ + yl = I. Finally, if we define fcg = |rg|“^ and ki = |ri|“^, we 
obtain the following two restrictions: 

(21) qx := + lizSxy — A:g = 0, 

(22) q2 := + liz{xiSiy - yiSix) - fci = 0. 

From Eq. the velocities and wi are 

(23) Wq — Si -f X rg — (si 3 ,, S\y -\- l\z^ 0)j 

(24) Wl = Si + li X Ti = (sij, — Sij, + 0), 

and the impulses Ag and Ai are given by 


(25) 

(26) 


Ag = Iwg - Wo\ = Ksix - WOxySly + hz “ W0y,-1C0z)|, 
Ai = lirj - Wl I = I {wl^ - S 12 , + lizyi , w* - si^ - hzXi 


so the cost function g = /2 = A§ + Af i 


IS 


q ■— (sia: Wqx) ~h (^ly “t“ Wz ^Oy) 4“ (iCgz) ~\~ (Wi,j. Six 4“ ^Iz^l) 4“ {^ly ^ly ^IzXi') 4“ {^'^iz'} ■ 

We compute the critical points of q using Eq. ( p0| ). In this case, V* =% because 

3(gi,g2) 


d{six,siy) 


~ WzVi ~ 0 
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is impossible since hz ^ 0 and y\ ^ 0. Therefore, Eq. ( |20[ ) reduces to 
(27) 


\rcrit _ 


qi = q2 = 


d{q,qi,q2) 




= 0 


We will solve Eqs. (271 using the technique explained in O Ch.2]. In order to do that, we com¬ 
puted the Grdbner basis of in the polynomial ring K[six,siy,liz] over the field of fractions 

K = Fiac {Q[ko,ki,xi,yi,wo,wl]/{xi + yf — 1)) with respect to the lexicographic monomial order 
Six > Sly > hz, obtaining I = (pi,P 2 ,P 3 ), where 

Pi = koyi ■ Six - (koXi - ki) ■ siy - (fco - ki) ■ hz 
P2 = (2(fco - klf + Sklkfyl) ■ Siy 

+ (d/cpTi -|- 2fcg?/2 ~ 4 A:q -I- ikQkiXiyf — Sk^kiXi — Sk^kiyl + Skgki 
+ ikokfxi + 2kQklyl - dfeg/ci) ■ llz 

+ {klxiyfwiy - k^xiyiivox - k^xiyuvl^ - klyfwlx - ^o^i^^iy + k^yiivox 

+ koyiwl^ - 2klkiXiyfwl^ - 2klkiXiylwly -b 2klkiXiyiWQx + ‘^■klkiXiyiwl^ 

- 2klkiylwly + 2klkiylwQx + 2klkiylwlx + ‘^klkiyl'wly - 2klkiyiWQx 

- ‘^klkiyiwl^ + k^klxiylwly - koklxiyiWQx - koklxiyiwl^ - koklylw\^ 

- koklylwly -b koklyiwox + koklyiwlx) ■ ijz 


-b (2 A:q -b Sk^kjyf - Ak^k^ + 2k{) ■ hz 


kt) 

u4„,2„ 


- (fcoTij/iicox + k^xiyiwlx + k^ylivoy -b k^yj-wly + 2kQkiXiylwoy + 2k^kiXiylwly 

- ‘^klkiylwQx - ‘^klkiylwl^ + fcp/cij/iWox + klkiyiw\^ - k^k^xiyiivox 

- klklxiyiwlx + klklylwoy + klklylwly - koklyiWox - kokfyiwU 
P3 = 21/1 • itz + (xiytKy - Xivlwox + Xiyfw*^ - ylwh + yfwly - yfwox + 

- (fcoa;ij/iii;ox + koXiyfwlx - 2koXiylwoy - 2koXiylwly - 2koXiyiWox 

- 2koXiyiwh + koyfwoy + koyfwly + 2koyfwox + 2koyfwl^ - 2koylwoy 

- 2kQylwly - 2koyiWox - 2fco?/iU’L -b 2fciaiij/iicox + 2kiXiyiwlx: - kiyfwox 

- kiyfwlx + 2kiyiWox + 2kiyiwlx) ■ hz 

- {Ak^xi — 2fcg?/2 + d/cg -b AkokiXiyf — SkokiXi + Skokiyl — Skoki + Ak^Xi — 2kly\ + Ak\) 

The equation pg allows one to solve for hz, which can be substituted in p 2 to get siy and, finally, in pi 
to obtain Six- This can be done since the leading coefficient in each equation is not zero. 

Now we deal with the case when the vectors rg and rq are linearly dependent. We can reduce to either 
To = 7^1 = (1: 0,0) or rg = —ri = (1, 0, 0). There is no need to use all the machinery that we developed 
so far to handle these two degenerate cases. The following discussion shows how to solve both situations 
with simple geometric arguments. 


In the former case, i.e. To = G = (1,0,0), we must have fcg = ki and iCg = Wi by Eqs. (2I)-(24|. 
The cost function /2 can be expressed entirely in terms of the independent variables w^y, Wg^, as 
follows: 


/2 = (wox - + {woy - Wg f + {woz - + (w*^ - + {wh - Wo f -b {w^ “ 
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The critical points can be found by setting the partial derivatives of /2 with respect to Wq^, w^y, Wq^ to 
zero and solving the resulting system of equations. Doing so, only one solution appears: 


Wox = Wlx = 


WOx + wL 


■; '^Oy — 


WOy + Wly 


, Wqz = Wlz = 


WOz + WL 


In the other case, i.e. ro = —fq = (1,0,0), the values of kg and ki are not necessarily equal, hence 
fo = (kg^, 0,0) and Ti = (—k^^, 0, 0), but wi can be expressed in terms of Wg. Indeed, the conservation 
laws for the angular momentum h and eccentricity vector e in the intermediate orbit imply that foXWg = 
fi X wi and Wg x (fo x Wg) — rg = wi x (fq x wi) — rq, respectively, from which it follows that: 

* k^ ,1, kj^ J, 

Wlx=Wg^, Wly = -—Wgy, Wl^=-—Wg^. 

Kg Kq 

The cost function /2 can be written in terms of the independent variables Wg^, Wgy, Wg^, as follows: 


ki 


/2 = (WOx - W*g^y + (WOy - W*gyy + (wo^ - W*gj'^ + (wj,, - Wg^f + + ^^gy^ + + 

Taking partial derivatives and solving the resulting system of equations, we get the unique solution: 


^Ox = 


WOx +Wl^ * _ 1 ^OWOy ~ ♦ _ , kgWg^ - fciW*,, 

^Oy ^0 ,_2 , ,_9 ) '^Oz ^0 


kg + ki 


kg + ki 


5. Optimality of the Hohmann transfer 

In this problem, we want to find the optimal 2-impulse transfer between concentric and coplanar circular 
orbits. Assuming that the plane that contains both initial and final orbits is orthogonal to (0, 0,1) and 
that the initial point is on the x-axis, we can reduce to the following situation: 


^0 — (0: 0; ^Oz)) h — (0, Oj ^2z)) ^0 — ^2 — (0, 0, 0), Tg — (1, 0, 0), 

where Iqz and l 2 z are not zero. 

The nine unknowns are the components of the vectors /i = (lix,liy,liz)^ si = (sia,, siy, siz) and iq = 
(xi, yi,Zi). The seven equations relating them are: 

' + WySly + WzSlz = 0 

Wx = 0 

+ hyVl + hzZl = 0 

< hzZi = 0 

2,2,2 1 
+ J/i + fo = 1 

^Oz = ^Ix + ^ly + ^Iz + Sly^lz “ Slz^ly 

, ^2z — ^Ix ^ly fo ^Iz fo ^l('^ly^lz '^’Iz^ly) fo yi(^lz^la: -^Ix^lz) fo Zi(six^ly 

Since l 2 z 7 ^ 0, we have that Zi = 0. Substituting lix = 0 and = 0 in the third equation, we get 
hyVi = 0. We will discuss two cases: liy = 0 and liy ^ 0 (which means that j/i = 0). 


Case liy = 0. Here we can reduce to li = (0,0, liz), si = (sixfoiyjO) and rq = (xi,j/i,0), subject to the 
following conditions: 

(x1 + y1 = l 
< ^Iz = ^Iz fo Sly^lz 

V ^ 2 z — ^Iz fo ^l^ly^lz yi’^lxWz 
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Since we want to minimize the function /i, we need to introduce the two extra variables Aq and Ai, 
which represent the magnitude of each impulse and the following two restrictions: 

f Ag = sJa; + (®ly + ~ ^Oz)^ 

1 A^ = (sia^ — lizUi + hzVl)^ + {Sly + hzXl — hzXl)^ 


Using Lagrange multipliers as in Eq.(19l, we computed the reduced Grobner basis in the polynomial 
ring Q( 1 o 2 ,/ 2 z)[Ai, ..., A 5 , Ag, Ai, sia;, Sly, Uz, a:i, i/i] with respect to the lexicographic monomial order 
Ai > • • • > A 5 > Ai > Ag > Six > Sly > liz > xi > 2 / 1 , obtaining the following two solutions: 


xi = - 1 , yi = 0 , hz = ± 


IL + 


2z 


Six — 0, Sly — 


72 _ 72 

‘Oz Uz J 

12 , 72 ‘1 

‘Oz > Uz 


The cost function /i = Ag + Ai at each of those solutions becomes: 


(28) 


/i = 




\/^0z + ^2 


-I 


Oz 


±V2ll 


^/F 


Oz 


P2z 


-I 


2z 


A simple computation shows that the sign of the optimum liz (and also the sign of the numerators in the 
previous expression) coincides with the sign of Iqz + hz- This case corresponds to the classical Hohmann 
solution. 


Case liy 7 ^ 0. Since hyyi = 0 and x^ + yl = 1, this case is only possible when yi = 0 and Xi = ±1. 
When Xi = 1, we have Ig^ = l^z, so the the orbits are of the same radius. If Iqz = hz, the initial and final 
orbits are exactly the same and no maneuver is needed. On the other hand, if Iqz = —hz, the satellite 
must change the direction of rotation in two impulses, both at the same point, so rg = ry and Wq = wi. 
There are clearly infinitely many optimal solutions with fi = Ai + A 2 = 2 | 7 i;g|. 

It only remains the case xi = —1. Here the unknowns are h = {0,liy,liz) and si = (sia;, siy, S 12 ), 
subject to the equations 

hySly + hzSlz = 0 


^ 1 ... “f ^i2 T Siyhz Sizhy 


'■Oz 

;2 


ly 

72 


hz — hy ^Iz Siyhz Sizl 


ly 


The impulses are 


^0 “ ^ix + i^iy + hz — hzY + (siz ~ hyY 

^1 = sfx + {Sly — hz + hzY + (siz + ^ly)^ 


and the cost function is fi = Ag + Ai. We computed the Grobner basis of Eq.(19l in the ring 


, ^ 2 z)[Ai, ..., A 5 , Ag, Ai, Six, Sly, siz, hy, hz] with respect to the lexicographic monomial order Ai > 
> A 5 > Ai > Ai > Six > Sly > Siz > hy > hz, obtaining the following two optimal solutions: 


, ^Oz + ^Ozhz 

hz — ~—“ 




■ 4^L^2z + hzhz + 


^hzhzhoz + hzhz + hz) 


hy = ±' 


jllz+llz 


/ 2 

hz 


Siz — 


72 _ 72 

‘Oz hz 


hz 


72 

hz 


72 _ /2 
_ hz hz 1 
~ 72 I 72 hz 
‘'fir “T to 


hz 


Six = 0 
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These solutions are defined only when m < < 02 , where oi and 02 are the real roots of a^ + 2a^ + 2a + 

1 = 0. In particular, the sign of must be negative. Substituting these solutions in the cost function 
fi and comparing with Eq.(28), it can be checked that the solution of the previous case is always better. 

6. Two ROTATED ELLIPSES 


In this orbit-to-orbit transfer problem, we will assume that the initial and final orbits are two identical 
ellipses rotated an angle a G ( 0 , 7 r] lying on the same plane. We will restrict our optimization to 
intermediate orbits that also lie within the same plane, i.e. to a two-dimensional orbit transfer problem. 
Without loss of generality, we can assume that Iq = (0,0, h = (0,0,/ 2 z), Sq = (soa,, soy, 0), S 2 = 
(s 2 x,S 2 y, 0 ) are given, and that we have to find fo = {xo,yo,0), fi = {xi,yi,0), h = ( 0 , 0 , hz) and 
= (sia,. Sly, 0). In order to guarantee that the initial and final orbits have the same eccentricity and 
semi-major axis, and to maximize the symmetry of the equations, we impose I 2 Z = ^Oz = 1 (since the 
problem does not depend on the semi-major axis), S 2 x = —sox and S 2 y = soy. This ensures that the 
orbits are identical, but rotated an angle a = 2 arctan(sox/soy). Both orbits are also symmetric with 
respect to the x-axis. Finally, we need two extra variables Aq and Ai to represent the two impulses. 
For general ellipses, with arbitrary semi-latus rectump, all the values h, Si, Wq, wi, Aq and Ai calculated 
in this section have to be divided by 

The discussion above reduces the problem to two parameters sqx, Soy, nine unknowns Xq, yo, 2 : 1 , yi, Sia,, 
Sly, hz, Aq, Ai, six equations 

' eqi :=xl -b yo = 1 
eq 2 :=xl + yf = 1 

eys .— h,z “b ^lz( 2 ^ 0 ^ 1 y yo^lx) 1 2 JoSoy -b yo^Oa, — 0 
< 2 

ey 4 .— h,z “b hz{^l^ly yi-Sla,) 1 XiSOy yiSOa, — 0 

eys - = ^0 = ('®0a: ~ Six)"^ + (SQy — Siy)^ -b (1 — hz)'^ + 2(1 — hz){xo{soy — Sly) — yo(sOx ~ Six)) 

^ eye :=Aj = (sqx + sia,)^ -b (soy — siy)^ -b (1 — hz)'^ + 2(1 — ^iz)(xi(soy — siy) -b yi(sox + six)) 

and a cost function fi = Aq -b Ai. 

After introducing the Lagrange multipliers, the algebraic problem has 15 equations and 15 unknowns. 
Although such a system is expected to have a finite number of solutions, this is not true in our problem, 
so some special treatment is needed. We will divide the problem in several cases, which will be discussed 
below. 


Case 1: We impose the extra condition yo -b yi ^ 0, which is done algebraically by introducing an 
additional variable k and adding the equation 1 — fc(yo -b yi) to the system. Geometrically, this new 
system looks for orbit transfers that are not symmetric with respect to the x-axis. We have no proof that 
the system has always a finite number of solutions, but we have collected extensive numerical evidence 
that this is indeed true. The best orbit transfer never happened to come from this case, as shown in 
Subsection lfi.il 

Case 2: Now we consider the remaining case, i.e. yo +yi = 0. It follows from eyi and ey 2 that xi = ±xo, 
so we split the analysis again: xq = xi (case 2a) and xq = —xi (case 2b). The former represents transfers 
whose initial and final points are symmetric with respect the x-axis, and the latter is a degenerate case 
when the initial point, the final point and the origin are collinear. Both cases have a finite number of 
solutions, which we will compute explicitly below. 

Case 2a: We assume here that Xi = Xq and yi = —yo- Subtracting ey 4 from eq^, we obtain that 
yoSix = 0. When yo = 0, we have Xq = Xi = ±1 and the cost function can be written as fi = 
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\/(sia: — + A + \/(six + sqxY + -^i where A is an expression that does not involve Since 

Six vanishes from all the equations when j/q = Oj we can consider it as a free variable. Setting the 
derivative of fi with respect to six to zero and solving the equation gives six = 0 after some algebraic 
manipulation. Now substituting xq = xi = ±1, yo = Ui = 0, six = 0 in the equations, leaves us 
with only one restriction ± hzSiy =F sqj, = 1. Solving for siy and substituting everything in the cost 
function fi, we get a minimization problem with only liz as a free variable, which gives the following 
two solutions: 


(29) Xq — Xi — dll, yo — yi — 0, Six — 0, Siy — SOy, llz — 1, fl — 2|soa:|- 

It only remains to see what happens when six = 0 and yo ^ 0. Comparing eq^ and eqo shows that 
Ao = Ai, so we can replace the cost function /i by 

2/2 = = ^Ox + i^Oy — Sly)^ + (1 ~ hz)'^ + 2(1 — ^lz)(a;o('SOy ~ Sly) ~ yoSOx)- 

This reduces the problem to four unknowns xq, yo, siyj hz, subject to two equations eqi and eqo- The 
case a;o = 0 leads easily to the following two solutions: 

^0 ~ ~ 9, yo — dll, yi = ^ix — O 5 siy = soy, liz — \^1 T ^ox, 

fl = 2|l=FS0a;T VlTSoxI- 


A straightforward verification shows that |soa:| > 2|l=psoa;T'\/l T soa;| for all sqx S (—1,1), which means 
that the solution (30) is always better than (291. 

From now on, we assume that xo ^ 0. This allows us to express siy in terms of xo, yo, hz using eqs, as 
follows: 

1 “t“ XoSoy yo^ox hz 


Sly - 


hzXo 


Substituting the expression for siy in the cost function \f2, we obtain a rational function c{xo,yo,hz)- 
Therefore, we have to minimize c subject to aig + yg — 1 = 0, which is equivalent to solving the equations 


dc 

'dyo 


dc 

Wz 

dc 

) o 

OXq 


= 0 
= 0 


a^o + 2/0 - 1 = 0 


We can clear denominators, without losing any information, by multiplying the first equation by ^f^Xg 
and the second one by l\zXo, since hz and Xq are both non-zero. We define the equations 

eg? = 0 

23/ dc dc \ 

=0 

SO our system is equivalent to solving eqi, eq-j and eq^. To solve these algebraic equations, we use 
resultants. Taking advantage of the fact that egi does not contain any term involving hz, we define 
P7,8 := R-eszi,(eg7,eg8) G 'L{sox, soy)[xo,yo\ and := Res3,o(egi,p7,8) G I.{sox, soy)[yo]- Any solution 

of our system satisfies both p 7_8 and Pi, 7 , 8 - Conversely, a solution of pi, 7,8 = 0 can be extended to a 
solution of the original system {egi, eg 7 , eg8} using the following procedure: 

(1) Find a solution yg G K of Pi, 7 , 8 ( 2 / 0 ) = 0. 
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(2) Write P7,8ixo,yo) = J2^>o^^iyo)xo G Z{sox, soy)[yo][xo] and substitute every even power 

by (1 — y^y and every odd power Xq^^^ by a;o(l — yoY- This way we obtain a polynomial of 
degree one in xq with coefficients in Z{sox, soy)[?/o]- This substitution is correct since Xq satisfies 
eqi = x^ + y^ - 1 = 0. 

(3) Assume that the polynomial obtained in the previous step is qi(yo)xo + qo(yo)- Then, calculate 
xo = -qo{yo)/qi{yo) e K. 

(4) Apply the Euclidean algorithm to the polynomials of eqy and eqg until a polynomial of degree 
one in liz appears. 

(5) Assume that the polynomial obtained in the previous step is ri{xQ,yo)liz + ro(xo,yo)- Then, 
calculate hz = -ro{xo,yo)/ri{xo,yo). 

The theory of resultants (see for instance Section 3 of 0) guarantees that the polynomials qi{yo) and 
ri(a;o,?/o) are different from zero. Since Pi, 7,8 is a polynomial of degree 48, there are potentially 48 
different solutions to our system of equations. A closer inspection of Piy.s shows that it factorizes as 

y^yo - l)®(2/o + y^ivlslx - 2yoSox + y^sly + 1 - sly^pohaiyo), 

where pol 2 oiyo) G Z(sox, soy)[j/o] is an irreducible polynomial of degree 20. Therefore, the roots of pi, 7,8 
are 0 with multiplicity 8, 1 and —1 with multiplicity 6, the complex numbers 


SO: 


SQy 


Ox 




- 1 


^Ox 




with multiplicity 4, and the 20 different roots of pol 2 o{yo) = 0. The roots 0, 1, —1 are discarded since we 
have already excluded these subcases. The complex roots with multiplicity 4 are not real since Sq^, + Sq^ 
is the square of the eccentricity, which is always < 1, by assumption. This leaves only the 20 roots of 
pohoivo) to be considered. 

Note also that the polynomials P 0 I 20 , qo, qi, ?'o and ri can be precomputed symbolically (as explained 
above), so the solution is: 


(31) 


2/0 = -2/1 = a root of po/20, xq = Xi = - 

1 “t“ XQSQy yO^Ox ^Iz 


Sly — 


WzXfl 


hz — — 


gi(2/o) 

90(2/0)’ 

T'i(a;o,2/o) 


Six — bi 


roixo,yo)' 


All together, we have 22 different solutions of case 2a: 2 from Eq. (30) and 20 from Eq. (31). Extensive 
numerical evidence shows that the best transfer always comes from one of these solutions, as discussed 
in Subsection HHl 


Case 2b: We assume here that yi = —yo and xi = —xq. In this case, we have \{eq8 + eqi) = 

1_^2 

— 1 + yiiSiix = 0, which implies that Wg = —la tli6 particular case when liz = 1, the following 
two solutions are found directly from the equations: 

(32) Xq — dll, Xi — ffil, 2/0 — 2/1 — b, llz — 1, Sly — SQyj Six G [ |soa,|, I-Soxl], fl — 2|soa,|. 


When llz = —1, then only one solution is possible: 


(33) 


^0 — 1; yo — yi — O7 ^ly — ^Oyi 

hz = -1, /l = (|1 + SOyl + |1 - SOyl) \J^ + S^x ■ 


^Ox^Oyi 


The solution (33) has fi > 2-\/4 + > 2|soa;|, so it is always worse that (32) and can be discarded. It 

can also be shown that solution (32) is always worse than (30), so it can be safely ignored as well. 
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The only remaining case is \liz\ ^ 1- Using eq^ — eg^ 4 , we can write six in terms of the other unknowns: 

'50y)'^0a: 


■Six — 


The problem has now only three variables xq, siy and hz and a single constraint eqi, which after 
1_;2 

substituting yo = becomes: 


eq,:= - 1) + {1 - = 0. 

The cost function fi = \/~S^ + \/^i is the sum of the square roots of two rational expressions in Xq, siy 
and liz- In this case, we will not introduce the extra variables Aq and Ai, since the square roots can be 
removed with an algebraic trick. First, according to the theory of Lagrange multipliers, we should have 
V/i = AVegg at each local extrema of /i. This produces the following equations: 

5A? 


(34) 


1 


5A2 


+ 


1 


2y^Qdsiy 2y^^dsiy 

1 9Ag 1 aAf 

2yA2 9xo 2y/Af dxQ 

1 9Ag 1 dl\l 

2y/^ dhz 2^^ dhz 


= 0 

“ 2Asg,jXo 

= -AXhzil-ll) 


We can remove A by multiplying the last two equations of (34) by 2 / 12(1 — and So 3 ;Xo, respectively, 
and adding them. 


(35) 


1 


2\/^ 


(9A^ 

2/12(1 - + SQx^O 


9A2\ 


1 


dhzj 2yAf 


9A^ 

2/lz(l — ^ ®0a;^0 


dAj 


= 0 


Finally, we remove the square roots in the first equation of Eq. (34) and also in (35) by moving one of the 
terms to the right, and then squaring both sides. After some algebraic manipulation and the introduction 
of the non-zero factors Ifzihz + l)^(^i 2 — 1)^ and /i 2 (/iz + l)(^i 2 ~ 1)^ to clear the denominators, we get 
the following two polynomials: 


6910 ^iz(^iz + ~ 1)^ 

6911 := ^izihz + l)(^i2 — 1)^ 


9A2 


ds 


1?/ 


A?- 


9Afy 

dSlyJ 


A^ 


= 0 




j aAgy 

^0x^0 -- I 


' £57 / 

ohz J 


Ai 


BAl 


( 2 / 12(1 / 12 ) + Sqj-Xq 


dhz) 


A^ 


= 0 


At this point we have reduced the whole case 2b to three polynomial equations {egg, egio, egn} in three 
unknowns xg, siy and / 12 . To solve the system, we exploit the fact that egg does not contain the 
variable Si^. Define pio,ii := ReSsj^ (egio, egn) G Z(soa:, Soi/)[a;o, ^iz] and pg, 10,11 := Res 2 ;o(egg,pio,ii) G 

) [^Iz] ■ 

A similar procedure to the one used in case 2a allows us to obtain a full solution of the equations 
{egg, egio, egn} from a zero of pg,io,ii- The procedure is described below: 

(1) Calculate a root hz G M of pg,io,ii = 0. Although pg,io,ii is a polynomial of degree 166 in / 12 , 
it can be factored as a product of polynomials that are either not zero by assumption or that 
have degree lower than 6. 
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(2) Use egg to obtain the two possible values of xg: 


Xo = ±\ 1 - 


1-^L 


If these values are not real numbers, then the following steps can be skipped. 

(3) Apply the Euclidean algorithm to the polynomials of egio and egn to compute their greatest 
common divisor. The algorithm stops when a polynomial of degree 1 in siy appears. 

(4) Assume that the polynomial obtained in the previous step is Ui{xQ,liz)siy + uq{xo,Iiz)- Then 
calculate si„ = 

The value of the remaining variables are: 


( 5 ) 


xi = -xo, yo = -yi = 


1 - /2 


— 




SQx Wz) 

The polynomials P9,io,ii’ '^o ui can be precomputed symbolically, following the same procedure as 
above. These polynomials can therefore be reused to calculate the solutions of case 2b for any given soa: 
and soy. 


(36) 


hz = a root of pg, 10 . 11 , xo = -xi = ±Wl - 


1 - 


, yo = -yi = 


1 - /2 
nz 

^Ox 


Sly — 


Uo{xo,liz) 


Wl(xo,Uz)’ 

We have collected extensive numerical evidence showing that these solutions are always worse than those 
of case 2a. Anyways, since there are only a finite number of solutions in this case, which can be computed 


XQ{llz^ly S0y)s0x 


by the explicit formula Eq. (36), we recommend that these solutions are included when looking for the 
best orbit transfer. 


6.1. Numerical tests. In our numerical computations, we explored a wide range of values of Sqx 
S oy, in such a way that all possible eccentricities and angles between the ellipses were considered. 


and 


When the ellipses are rotated 180 degrees, the optimal solution is always provided by Eq. (30) of case2a. 
Indeed, easel does not have a real solution if the eccentricity is 0.1,0.2, 
given by Eq. (36) and case2a given by Eq. (31) exist but are worse. 


,0.9. Solutions of case2b 


In the rest of the cases, the solution of case2a given by Eq. (31) is always the best. To check this 
efficiently, we used only rational values for sqx and soy The trick to achieve this is to set 

2ab 


^Ox — 6 


+ &2’ 


^Oy — ^ 


62 


where e is the desired eccentricity and a and b are integers chosen in such a way that the desired angle 

( 2 _7 2 


2ab 


In the following tests, we used e = 0.1,0.2,..., 0.9 and the pairs (a, b) were selected to approximate the 
angles a = 5,10,..., 175 degrees. The case a = 180 was discussed above. 

For each value of a, b and e, we computed the best solution of case 1 (solving the system numerically). 


case2a using both Eq. (31) and Eq. (30), and case2b using Eq. (36). In total, we have explored more 


than a thousand test cases. We extracted several conclusions from the data we computed. 

First of all, the solution of case2a is indeed the best one, as we mentioned before. In Figure we show 
the fuel of this transfer (or rather, the cost function /i = Aq + Ai, which is proportional to it) obtained 
as a function of e and a, and what this transfer would look like when e = 0.7 and the angle between the 
ellipses is 85 degrees. 
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DeltaJ + Delta_1 of best transfer Besf orbif transfer for ecc =0.7 and angle =85 deg 




Figure 2. Aq + Ai of the best transfer and an example of one of these transfers. 

We can also compare the angle between the semi-major axis of the initial orbit and the direction given 
by tq. We observe in Figure [^a) that when a is small (less than 40 degrees), this separation is higher 
than 50 degrees for eccentricities up to 0.5. When a is near 180 degrees, this separation becomes smaller 
(and is zero in the case a = 180). 

The separation shown above led us to study how much fuel can be saved by using our optimal transfer 
instead of the one from apogee to apogee. Figure [^b) shows the ratio (in percentage) between the fuel 
consumption of both transfers. 


Angle between semi-major axis of initial orbit and (x_0,y_0) 



Fuel comparison: Best transfer vs Best apogee-to-apogee transfer 



(a) 


(b) 


Figure 3. (a) Angle of separation between the semi-major axis of the initial orbit and 
the point where the first impulse is applied, (b) Fuel comparison between the best 
transfer and the best one from apogee to apogee. 


On the other hand, easel does not always produce a valid real solution. Even in those situations where 
easel provides a solution, it is always very poor compared to the ones of case2a, up to one or two orders 
of magnitude worse depending on the eccentricity. Case2b always produces valid solutions, but they are 
as bad as those of easel. 


Finally, the solutions of case2a given by Eq. (30) are worse than the optimal one, but they are no more 


than 10% worse for eccentricities below 0.6 and up to a 55% worse for higher eccentricities, as shown in 
Figure 
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Fuel comparison: Best transfer vs Solutions ot Eq.(30) 



Figure 4. Fuel comparison between the best transfer and the solutions provided by Eq.(30) 


7. Conclusions 

Firstly, in Section]^ we have presented an innovative approach to the study of the multi-impulse orbit 
transfer problem with minimum fuel, that only requires to solve a system of polynomial equations. No 
trigonometrical functions are needed. This allows one to use all the algebraic machinery that is currently 
available: Grobner bases. Resultants, Elimination theory. Root Isolation, etc., which was not possible 
with the previous formulations of the problem. 

In Section we provided an alternative method for solving the point-to-point two-dimensional orbit 
transfer problem studied in [T] by Avendaho and Mortari. The significant advantage of our new approach 
is the efficiency of the computation, since only arithmetic operations are used. 

Moreover, in Section]^ we used the well-known Hohmann transfer problem between two circular orbits, 
to show the power of our technique. We worked under very general assumptions, i.e. letting the transfer 
orbit to be out of plane, but we showed that the best transfer is indeed coplanar. The novelty of our 
analysis is that we allowed the initial and final orbits to have angular momentum pointing in opposite 
directions. Even for such an extreme case, the classical solution is proven optimal. 

Furthermore, in Sectionwe analyzed in depth the problem of changing between two identical elliptical 
orbits of eccentricity e which are coplanar and rotated a certain angle a. Here we restricted our search 
to transfer orbits that are coplanar with the other two. The first surprising result that we got is that 
the optimal transfer does not go from apogee to apogee (except when a = 180deg), but it is separated 
from the apogee a certain variable angle that depends on e and a. This angle is higher than 50 degrees 
when e < 0.5 and a < 40 deg. For lower eccentricities and small values of a, the separation angle can 
be as large as 80 degrees. 

The large difference between the best transfer orbit and the best one from apogee to apogee means also 
a signihcant difference in the fuel consumption of both maneuvers. For angles a up to 80 degrees, the 
savings obtained by using our transfer orbit are always higher than 25%, and for angles a < 10 deg, our 
transfer consumes less than half of the fuel needed to go from apogee to apogee. 

Finally, solving the problem has been reduced to the study of several subproblems, each of which consists 
of a set of polynomial equations. All but one of the subproblems have been solved symbolically, which 
means that an explicit solution is available given any initial and final orbits. The remaining subproblem 
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can be solved numerically for any given data, but the optimal solution does not seem to come from this 
particular subcase. Indeed, numerical evidence suggests that the optimal solution always comes from 
the same subcase, for which we have a symbolic solution. 
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